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The superfluid-insulator transition in the disordered two-dimensional Bose-Hubbard 
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We investigate the superfluid-insulator transition in the disordered two-dimensional Bose-Hubbard 
model through quantum Monte Carlo simulations. The Bose-Hubbard model is studied in the 
presence of site disorder and the quantum critical point between the Bose-glass and superfluid is 
determined in the grand canonical ensemble at fi/U = (close to p = 0.5), n/U — 0.375 (close to 
p = 1), and fi/U = 1 as well as in the canonical ensemble at p = 0.5 and 1. Particular attention is 
paid to disorder averaging and it is shown that a large number of disorder realizations are needed in 
order to obtain reliable results. Typically, more than 100, 000 disorder realizations were used. In the 
grand canonical ensemble, we find Ztc/U — 0.112(1) with n/U = 0.375, significantly different from 
previous studies. When compared to the critical point in the absence of disorder {Ztc/U = 0.2385), 
this result confirms previous findings showing that disorder enlarges the superfluid region. At the 
critical point, we then study the dynamic conductivity. 

PACS numbers: 61.43.Bn, 05.60.Gg, 05.70.Jk, 02.70.Ss 



I. INTRODUCTION 



Since the publication of seminal papers on superfluid- 
insulator transition by Fisher et al. ,^''^ the disordered 
two-dimensional (2D) Bose-Hubbard (BH) model has at- 
tracted much theoretical attention. The disordered BH 
model with Cooper pairs acting as charge-2e bosons has 
been argued to describe the superconductor-insulator 
transition in thin amorphous films3""^ Recently, with the 
development of experimental techniques for constructing 
the BH model by confining cold alkali atoms in optical 
lattices, superfluid (SF) to Mott insulator (MI) transi- 
tions have been observed in clean optical lattices^. Fur- 
ther experiments have introduced disorder in the optical 
lattices by speckle fields to investigate the phase diagram 
of a disordered three dimensional optical lattice^i^. Al- 
though experimental techniques still need to be refined, 
there is renewed interest in the disordered BH model from 
experiment. From theory, there is a long-standing re- 
search interest in studying the disorder induced phase 
transitions and searching for Bose glass (BG) phasei^i"— 
In particular, recent work has focused on the transition 
in 3 dimensions^li^"— and the validity of the relation 
z = d established by Fisher et alJ» has been questioned 
both in numericalMi^ and theoretical studies'^^ . Here we 
shall focus exclusively on the disordered two-dimensional 
(2D) model where a number of outstanding questions re- 
main. 
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FIG. 1: (Color online) The phase diagram of the 2D Bose- 
Hubbard model. Shaded areas are the Mott-insulating phases 
for zero disorder as determined from strong coupling expan- 
sions in Ref. [S^ls^. The mean field phase boundaries and 
constant density profiles for zero disorder are shown as red 
dotted lines. The dashed line indicates the constant chemi- 
cal potential fi/U = 0.375. The solid triangle indicates the 
location of the transition to the Mott phase in the absence of 
disorder as determined by SSE simulations along the dashed 
line from Ref 39. The three solid squares from bottom up 
are for the locations of superfluid to Bose glass transitions in 
the presence of disorder at n/U — 0, 0.375, 1, respectively, as 
determined from SSE simulations in the present work. 



The disordered 2D BH Hamiltonian is given by 



U 



H = -i^(at^^a,+H.c.) + -^n,(n,;~l)-|-^(e,-//) 
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(1) 

where (5 = x, y, i is the nearest neighbor hopping am- 
plitude, al{ai) is a boson creation (annihilation) opera- 



tor, H.c. means Hermitian conjugate, n^ = aja^ is the 
number operator, U is on-site interaction, fi is chemical 
potential and e^ is uniformly distributed in the interval 
[—A, A], with A controlling the disorder strength. Here- 
after, we shall explicitly give energy values for both t and 
U in each calculation to facilitate comparisons with cal- 
culations in the literature, since Hamiltonian definitions 
in each paper may be different. The phase diagram for 
this model in the absence of disorder as determined from 
strong coupling expansions (from Ref. I37ll38f ) is shown in 
Fig. [1] with Mott insulating lobes with fixed on-site par- 
ticle number extending into the superfluid phase. The 
mean field phase diagram^ is shown as dotted lines in 
Fig. [1] and the importance of fluctuations is clearly ev- 
ident from the discrepancy between the mean field and 
strong coupling results at the tip of the envelopes. De- 
spite its relative simplicity, a detailed understanding of 
this model with disorder has proven surprisingly difficult 
in particular for numerical work. 

Existing studies using various methods address differ- 
ent aspects of the disordered BH model, and often arrive 
at contradicting conclusions. It is therefore most useful 
to revisit this problem using current high performance 
numerical techniques. In the presence of disorder it is 
knowni that a BG phase appears in addition to the SF 
and MI phase present without disorder. The question of 
whether a transition directly from the SF to the MI with- 
out an intermediate BG phase is possible in the presence 
of disorder arises. However, this question now seems set- 
tled with a proof that there is always an intermediate BG 
phasci^ In the simulations we report here, we are always 
in the strong-disorder regime (A = U) and we focus on 
the SF-BG transition since we expect the MI phases to 
be strongly suppressed at strong disorder. 

A model closely related to the Bose-Hubbard model 
is the {N — 2) quantum rotors model: [cos(6'r), sin(0r)], 
believed to be in the same universality class as Eq. ([T|). 
This model describes a wide range of phase transitions 
dominated by phase-fiuctuations: 
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Here, t is the renormalized hopping strength and i^- = 
Ly. is the angular momentum of the quantum rotor. 
The angular momentum can be thought of as describ- 
ing the deviation of the particle number from its mean, 
Lr ~ rir — ng. This model can be obtained from Eq. ([T]) 
if amplitude fiuctuations are integrated out when com- 
pared to Eq. ([1]). It is implicitly assumed that only phase- 
fluctuations are important at the quantum critical point 
(QCP), while amplitude fiuctuations are neglected. The 
model Eq. ^ can be simulated very efficiently using the 
Villain (link-current) representation, and a specialized di- 
rected geometrical worm algorithm^^i?=§. has been devel- 
oped for this purpose. This technique can be applied to 
both the clean and disordered model and high-precision 
results can be obtained. However, when compared to 



results obtained by direct simulations of Eq. ([T]) by Ba- 
trouni et al. — as well as from simulations of Eq. ^ in 
the hard-core limit by Makivic et al. '*'*, discrepancies 
appeared in particular for the universal features of the 
conductivity at the critical point, as we discuss below. 
Here we shall show that a likely explanation for these 
discrepancies is an inaccurate determination of the QCP 
in the direct simulations of Eq. ([T]). Furthermore, it is 
important, from our point of view, to provide accurate 
quantum critical parameter values (e.g. tc and Uc for the 
Hamiltonian defined in Eq. ([Ij) using recently improved 
simulation techniques. 

Our calculations proceed under two different condi- 
tions: fixed particle number (canonical) in order to fol- 
low the calculations of Ref. '41' or fixed chemical potential 
(grand canonical) . This corresponds to taking two differ- 
ent routes in probing QCP in the BH model phase dia- 
gram illustrated by the dotted and dashed lines in Fig. [T] 
Previous calculations on the clean system by Smakov et 
al^ using a quantum Monte Carlo (QMC) method have 
located the QCP at Ztc/U = 0.2385 (solid triangle in 
Fig. [T]) at ^/U = 0.375 (equivalent io p — 1), where 
Z = 4 is the coordination number for the square lat- 
tice. As evident from Fig. [T] this is in excellent agree- 
ment with the strong-coupling result from Ref. I37ll38l . 
Here we follow a similar approach to determine the QCP 
but with an added on-site disordered chemical potentials 
Ei S [—A, A], with A = [/ in the region of strong disor- 
der. 

At the critical point in a 2D system, the dc conductiv- 
ity was predicted^ to have a universal value, a* , close 
to the conductivity "quantum" uq = e*"^ /h, with e* 
the charge of the bosons. However, the exact univer- 
sal value has yet to be determined. Experiments^!^ sug- 
gest the dc conductivity value to be a little larger than 
Uq, but there are concerns that the experimental tem- 
perature (typically > 0.5 K) is not low enough. Us- 
ing the quantum rotor model, Eq. ([2]), S0rensen et ali^ 
showed that the universal conductivity value is given by 
(T* = (0.14 ± 0.01)crQ. Previously, exact diagonalization 
on a hard-core BH model by Runge^ gave a value of 
a* = (0.15±0.01)cr(3 in good agreement. Later, by doing 
world-line quantum Monte Carlo (WLQMC) simulations 
directly on a disordered 2D BH model, Batrouni et al^ 
found a* = (0.45 ± 0.07)crQ and Makivic et al. ^ found 
a* = (1.2 ± 0.2)crQ using a hard-core BH model. These 
calculations differ from each other and the experimental 
value. It was later pointed ouf*"^ that in d spatial dimen- 
sions, the dynamic conductivity a{uj) obeys the following 
scaling relation near the QCP 



a(.) = 2wg(M:)(-^)AE(^), 
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which might in part explain the difference between the 
experimental value and numerical calculations since the 
fvjj /ksT — ^ is achieved in the experiments while in the 
imaginary-time QMC simulations with Matsubara fre- 
quency Wfc one always has huj^/kBT = 27rfc > 1. The 



regimes lThj/{kBT) <Si 1 and hcd/{kBT) ^ 1 are domi- 
nated by different transport mechanisms, hydrodynamic- 
collision dominated and collisionless phase-coherent, re- 
spectively. Hence, there is httle reason to beheve that 
a simple extrapolation using only huj /ksT ^ 1 can cor- 
rectly determine the observed experimental dc conduc- 
tivity. If a careful extrapolation first to L — >■ oo and then 
T — >■ is performed it is possible to gain some informa- 
tion about this limit <^ Nevertheless, the entire function 
S in Eq. ([3|) is universal and the above mentioned numer- 
ical results should still agree on this universal function 
and the extrapolations should yield the same number. 
However, that extrapolated number which we shall call 
a* might not be closely related to the dc conductivity 
but will instead correspond to a higher frequency part of 
S. A sketch of the expected behavior is shown in Fig. [21 
Our results here show that if a careful determination of 
the QCP is performed, Eq. ([T]) and Eq. ^ yield the same 
value for a*. 
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FIG. 2: A sketch of the conductivity with the two regimes, 
fko/{kBT) ^ 1 and hoj/{kBT) <C 1 clearly apparent. Also 
shown is the difference between the universal dc conductivity 
a* and the conductivity a* obtained from extrapolating the 
high-frequency numerical data. An approximate position of 
the first non-zero Matsubara frequency is also shown. 



If the conductivity for the BH model at the critical 
point is a scaling function of fiw/fc^T, one would expect 
plots of cr(w) versus uj to show deviations from scaling 
at the critical point even for finite systems. A primary 
goal of this paper is to see if such deviations are observ- 
able for the two-dimensional disordered BH model if the 
QCP is determined carefully. We believe that the results 
we present here show clear indications of such deviations 
from scaling with w. 

Not only the conductivity value and dynamic conduc- 
tivity scaling differ in historical studies, but also the QCP 
has diverse estimates. For example, Zhang et alM used 
both ground-state and finite-temperature QMC simula- 
tions to locate the QCP of a hard-core BH model, which 
seems to agree with earlier work of Krauth, Trivedi, and 



CeperleyiiS, These two simulations are, however, quite dif- 
ferent from what Batrouni et al. found^. It is notewor- 
thy that these QMC calculations have used only about 
100 disorder samples; we show below that this is not suf- 
ficient for a precise determination of the critical point. 

Recently, QMC algorithms for the BH model have un- 
dergone a big improvement in efficiency with the devel- 
opment of stochastic series expansion (SSE) algorithm^ 
and the directed loop-update technique4^i^ This algo- 
rithm does not have the imaginary-time discretization 
error inherent in the traditional WLQMC method. The 
loop updates are especially important in the quantum 
critical region where the long correlation time substan- 
tially increases the errors of the WLQMC simulations. 
It is now feasible to check the convergence of measured 
quantities with the disorder averaging by increasing the 
number of disorder samples dramatically. This has been 
found by Wallin et al.^^ and Hitchcock et al."^"^ to be cru- 
cial in order to obtain reliable results. 

In this paper, we first discuss the SSE algorithm as it 
is applied to the BH model. In particular, we discuss 
the measurement of superffuid density and dynamic con- 
ductivity. Then, tests on the equilibration process and 
autocorrelation functions of the simulations are described 
to validate the estimates of the superffuid density needed 
to locate the QCP. We then show the superffuid density 
scaling figures for both the canonical and grand canonical 
ensembles, discussing the connection and differences from 
previous results. Finally, we show dynamical conductiv- 
ity scaling and the extrapolated universal conductivity 
values in the high-frequency limit. 



II. NUMERICAL METHOD AND 
CONVERGENCE TESTS 

As pointed out by Weichman,^^ analytic solutions of 
the model, Eq. ([T]), based on perturbation of the non- 
interacting limit has a fundamental difficulty because in 
the absence of repulsive interaction U or the Pauli exclu- 
sion for fermions, the presence of any disorder no matter 
how weak will condense a macroscopic number of parti- 
cles into the lowest localized free particle eigenstate of 
the random potential. This is not a meaningful state to 
do perturbation on. When the interaction U is non-zero, 
a complicated competition between the interaction and 
disorder potential makes it impossible to do analytic cal- 
culations. In this paper, we resort to QMC simulations 
of the model with the SSE algorithm;^"—, since it is able 
to treat any interaction strength and disorder realization. 
The method is brieffy discussed below. 

In the SSE formalism, the above site representation 
of Hamiltonian Eq. ([T]) needs to be written in a bond 
representation. 



H 



b=l 



Ifc 



H2b) 



(4) 



where the b = (ij) is the bond index, Hn, is a diagonal 
operator, and i?26 is an off-diagonal operator given by 
the following equations: 



Him = C - [—ni{ni - 1) + —nj{nj - 1) + (e,; - p.)ni 



H2M = t{alaj +H.C.), 



(5) 



where A = A/Z, with A one of ^,ei,U and Z — 4 being 
the coordination number of each lattice site. iVf, — ZN/2 
is number of bonds in the system with N lattice sites. C 
is a constant chosen to ensure a positive definite expan- 
sion. Since e is a random number, we need to use the 
disorder amplitude A in order to determine an appropri- 
ate value for the constant C. 

The partition function can then be expanded as 



Z = Tr{e-^^} =J2i: 
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a n=0 S„ ' 4=1 

where {|a)} represents a complete state basis and 

Sn = [ai,&i], [a2,&2],---, [a„,6„] 

is an operator-bond sequence, with ai G {1,2} denoting 
the type of operator (l=diagonal, 2=ofi'-diagonal), and 
bi € {1, . . . , Nb} being the bond index. The expansion is 
now truncated at order M, and M — n unit operators for 
the nth order term are inserted: 



a Sa 



Here M has been made large enough so that the probabil- 
ity for the expansion order exceeding M can be neglected. 
The Monte Carlo weight for the configuration (a, Sm) is 
then given by: 



1^(0,3 m) 



M! ' 



The resulting configuration space, consisting of state and 
operator sequences, is updated with both diagonal and 
loop updates. 

The diagonal update probability with n ^ 71 ± 1 is 
given by 



P{n — > n - 1) = 



m(3{a{p)\HiM\a{p)) 

M -n 

M -n+1 



(9) 



mi3{aip)\HiM\aip)) 

where \a{p)) is the propagated state 

p 
Hp))^llHa„tAa). (10) 



For the loop update, transition probabilities during 
the loop construction are chosen to have bounce-free 
or bounce minimizing solutions of the operator vertex 
equation."^ The configuration weights for operator ver- 
tex in the disordered BH model are, for particle occu- 
pation Ui and rij for nearest neighbor (ij), given by the 
following different cases: 

w{ni,nj;ni + l,nj - 1) = Uj^Jirii + l)nj{l - Sm^n^^J 
w{ni,nj]n^ - l,nj + 1) = Uj^fiiin^ + 1)(1 - (5„3,n,„,x) 

w{ni,nj;ni,nj) = C - [— 
+ {ei - p,)ni + (ij - fi)nj]. 



w{ni,nj\ni,nj) = C - [—ni{ni - 1) + -rnj{nj - 1) 



(11) 



Here nmax is the maximum number of bosons that can 
occupy the same site, rimax is usually assigned a large 
enough value (compared with the average density of the 
system) so that it can allow particle number fluctuations 
while at the same time will never actually be exceeded 
during the QMC simulations. From the operator vertex 
weight expressions, we can see that for the disordered 
BH model, all these weights need to be calculated "on 
the fly". Note that we have not explored the possibil- 
ity of tabulating all the operator vertex weights, which 
could speedup the calculation. In the tabulation, one 
would need to set up ZN/2 probability tables (each ta- 
ble corresponds to one bond) ; the number of elements in 
the table is determined by rimax, e.g. for Timax — 4, each 
table has 3392 elements. 

The total energy of the system is related to the expan- 
sion order of the operator sequence 



E = 






(12) 



\{Ha^.t,M). (8) {W^) 



where n is number of non-unit operators in the operator 
sequence as discussed above. 

The physical quantity that indicates the superfluid- 
insulator transition is the normalized superfluid density 
Ps,*^ computed as the average square winding numbers 



2t(3p' 



(13) 



where W^ — W^ + Wy, (3 is inverse temperature, and p 
is the average number of particles per lattice site. 

To locate QCP, we will use the normalized superfluid 
density finite-size scaling relation^ 



L"/(aLi/''(5,/3i~^), 



(14) 



where L is the linear dimension of the square lattice, 
d = 2 is lattice dimension, a = 2 — d — z, 5 measures the 
distance to the critical point (e.g., to determine Uc, we 
have 5 — {U~Uc)/Uc, and similarly for i^), -^ is dynamical 
exponent, which is predicted^ to be z = 2, a is a non- 
universal metric number, and the function / is universal. 



Hence, we have a = — 2, and if we keep /3L~^ fixed, and 
plot LP'Ps versus Zt/U for different lattice sizes, all the 
curves will intersect at the critical value of Uc ■ Note that 
we assume the validity of the relation z = d which has 
been brought into question recently."^® 

The conductivity of the BH model as a function of 
Matsubura frequency oj^ = 2TTk/{hf3) can be calculated 
from the linear response relations (we illustrate using the 
X direction; similar response formulas apply for the y 
direction)^"^ 



aiiujk) = Sttctq 



(k^) - Axx{iu>k) 



Wfc 



(15) 



where {kx) is kinetic energy per link along the x direction, 
and A^x (ii-^k ) is the Fourier transform of the imaginary- 
time current-current correlation function Axx{t) 
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Here N is total number of lattice sites, and Axxij) = 
{jxiT)jx(fi)), with the paramagnetic current being given 

by jx{0) = itY.ri'^l+x'^r - alur+x), and Jx{t) being the 
Heisenberg representation of ja;(0). Note that the small- 
est Matsubara frequency wi corresponds to huji/{kBT) = 
27r significantly larger than one. Hence, as pointed out 
by Damle et al. ■^ this type of imaginary-time QMC cal- 
culation will necessarily be in the the coUisionless phase 
coherent regime with hw / [ksT) ^ 1. In the following 
we set fi = 1. 

To calculate the imaginary-time current-current corre- 
lation function Axx{t), we follow the discussion of Ref. 
Isg and divide Axx {t) into four components (which may 
be thought of as combinations of current and anti-current 
terms) by 



Axx{t) 
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A2^(r), 



(17) 



and 



A2^(r) = V(i^J(r,r)i^,-(O,0)), 



^^(i",t) = tot+^(r)ar(r), 
Kxi'^.T) = -tal{T)ar+x{T). 



(19) 
(20) 



Finally the imaginary-time current-current correlation 
function can be calculated with a binomial summation 
in QMC simulations^M^ 



Axx{i(^k) 
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7,^— ± 

n-2 

r;7,t/— ± m— 



A2^(t), 

_. n-2 

o E amn(iwfe)iV(zy,7;m)),(21) 



where amn{i<^k) is the degenerate hypergeometric (Rum- 
mer) function,^ i.e., Umnii^k) =iFi{m + l,n]il3ujk), 



N{v^ 7; m) is the number of times the operators KJ and 
K^ appear in the SSE operator sequence separated by ra 
operators, and n is the expansion order. The virtue of the 
above formulation is that the imaginary-time integral in 
Eq. (J16p is performed analytically, hence eliminating the 
discretization error of the imaginary-time integral that is 
present in the conventional WLQMC— 



1. Convergence Tests 




FIG. 3: The energy difference between two replicas Ea — E-y 
as a function of MCS averaged over 100 disorder samples. 
Results are shown for L — 16, t = 0.5, U — 11, and /3 — 12.8. 
We start recording the energy difference after 50 MCS. 
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FIG. 4: The autocorrelation function, C^Y2{ta), for W^ aver- 
aged over 100 disorder samples for L = 16, t = 0.5, U — 11, 
and /5 = 12.8. Note that the j/-axis is on a log scale. 



Numerically we perform SSE simulations for lattice 
sizes L = 6, 8, 10, 12, 16. According to scaling relation 
Eq. (|14p . the corresponding inverse temperatures for 



these lattices are /3 = 1.8,3.2,5.0,7.2,12.8, leading to 
a constant /3L~^ = 0.05. Our strategy is to ensure that 
each disorder realization has been equilibrated, resulting 
in a few statistically independent measurements for each 
disorder realization, with the error-bars then obtained 
from the disorder averaging where each realization can 
be considered statistically independent. We have per- 
formed several tests to ensure that this equilibration is 
attained and that the number of disorder realizations are 
sufficient21. 

We begin by considering the equilibration of the energy 
between two different replicas with the same disorder re- 
alization. We define one Monte Carlo sweep (MCS) as 1 
diagonal update, which includes a sweep through all diag- 
onal operators in the SSE expansion, followed by 10 loop 
updates (the number of loop updates included in 1 MCS 
is arbitrary, but it is usually determined by the ratio of 
the SSE expansion order to the average loop length). In 
Fig. [2J we show the equilibration of the energy difference 
between two energy replicas as a function of MCS for 
L = 16, t = 0.5, U = 11 averaged over 100 disorder real- 
izations. As mentioned, replica simulation means that for 
each disorder realization, we start two parallel SSE simu- 
lations a and 7 with vastly different initial configurations 
but with the same disorder realization, and monitor the 
evolution of Ea— E^. We see that after about 500 MCS, 
Ea — E~f fluctuates around zero, showing equilibration. 
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FIG. 5: (Color online) Convergence of the superfluid density 
for t = 0.5 and U = 11 as a function of the number of disorder 
samples A^sampic for various lattice sizes, L = 6, 8, 10, and 16. 

We then proceed check the convergence of the winding 
number W'^ by calculating its autocorrelation function. 



Cw^{ta) 
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(22) 



Our results for Cw^[ta) versus ta (MCS) are shown in 
Fig. |4]on a log scale averaged over 100 disorder samples. 
Since we expect the leading term in Cnf^{ta) to scale as 
, we estimate the autocorrelation time to be around 



-tjT 



■ 40 MCS. 



Based on the above tests we set the following parame- 
ters for the SSE simulations and measurements: for each 
disorder realization we perform 1000 MCS to warm up 
the system, which is followed by 1000 MCS measurements 
with each measurement being separated by 1 MCS. While 
this might seem inadequate to obtain small error bars 
for a single disorder realization, this approach is in fact 
optimal since each disorder realization is statistically in- 
dependent; reliable error-bars can then be obtained from 
the disorder averaging. We have done tests for the worst 
case, i.e., longest warm-up MCS and W"^ autocorrelation 
time for L = 16, t = 0.5, U — 11, and /? — 12.8, so we ex- 
pect the parameters to be sufficient for the other lattice 
settings. 

Finally, we focus on the disorder averaging which we 
denote by [•] to distinguish it from thermal averages 
< • >. Often one finds that the estimate of a quantity 
one wishes to calculate obeys a rather broad distribution 
with a substantial tail-- with a difference between the av- 
erage and the most likely values. This will be reflected in 
an underestimation of the error-bar calculated from the 
standard deviation and a slow "movement" of the av- 
erage as the number of disorder realizations is increased 
with the error-bar essentially constant with respect to the 
number of averages. We check our convergence here by 
focusing on the behavior of the superfluid density as we 
increase the number of disorder samples -/Vsampie- See 
Fig. [51 For each curve in the flgure, the first point 
is for 1,000 disorder realizations, the second for 10,000 
disorder realizations, and for subsequent points, the in- 
crement is 10,000 disorder samples. It is important to 
note that 1,000 disorder realizations is already an order 
of magnitude more than the previous calculations f^^^i^ 
where the number of disorder realizations was taken to 
be around 100. Here, we find that 1,000 disorder samples 
are still far from convergence. The idea of self-averaging 
does not seem to apply here, either. For the remainder 
of our results we have used in excess of 100,000 disorder 
realizations to ensure reliable results. 



III. RESULTS 
A. Fixed particle number, p — 0.5, 1 

As described, the SSE is essentially a grand canoni- 
cal simulation algorithm and suitable for fixed chemical 
potential QMC simulations. However, in order to be as 
close as possible to previous calculations by Batrouni et 
al. j**^ we also want to do simulations with fixed par- 
ticle number for each disorder realization. It is possi- 
ble to forbid particle number fluctuation during the loop 
construction, ^'^■^^ but we will adopt a simpler approach. 
We first do some rough estimates of chemical potentials 
to yield total particle number close to the required value 
for a set of disorder realizations. We then use the average 
chemical potential for all disorder realizations in the con- 
struction of the probability table during the loop update. 



and reject any closed loops that change the total parti- 
cle number. In practice, we find that with the roughly 
estimated chemical potential the loss of efiiciency caused 
by additional rejection is insignificant, typically 5% more 
rejection than the grand canonical simulations. 
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FIG. 6: (Color online) Superfluid density scaling from canon- 
ical simulations, i.e., fixed particle number with density 
p — 0.5, for A = 12 and t — 0.5, for various lattice sizes, 
L = 6,8,10, and 16. 



We now turn to our results using SSE simulations in 
the canonical ensemble with nmax — 4, fixed particle den- 
sity p — 0.5, and fixed disorder strength A = 12 and hop- 
ping integral t — 0.5, which are the same as in Ref. |4l|. 
Note that canonical systems with fixed particle density 
of p = 0.5 are very close to the grand canonical systems 
with p/U = 0. A series of U values will be used in the 
simulations to search for the critical point Uc based on 
finite-size scaling relation for the normalized superfluid 
density Eq. (fT4|) . 

To locate the QCP Uc or equivalently Zt/Uc, we plot in 
Fig. [5] the superfluid density L^ps as a function of Zt/U 
(fixing t = 0.5 and using a U grid of C/ = 11, 12, 13, and 
14) for various lattice sizes L = 6,8,10, and 16. All 4 
curves intersect very close to the point Zt/Uc — 0.182. 
The crossing of the curves is not as "perfect" as for the 
clean system where larger system sizes can be simulated 
and it would appear that corrections to scaling for the 
smallest size L = 6 might be sizable. However, for dis- 
ordered systems the results in Fig. [6] is at par with the 
best one can obtain. We also note that the fact that 
we do obtain a crossing of the curves at a single point 
is an indirect validation of the scaling relation Eq. p^ 
with z — d. Based on the data shown in Fig. |6l we es- 
timate the critical point to be Zt/Uc = 0.182 ± 0.003 
or C/c = 11.0 ± 0.5. To achieve a more accurate result 
will remain a challenge for some time. Since canonical 
simulations with p = 0.5 correspond closely to the grand 
canonical simulations with p/U = 0, we can compare the 
QCP obtained in these two simulations, and find that 
this canonical QCP value for p = 0.5 is very close to the 



grand canonical value of {Zt/U)c — 0.192 at fi/U — 0, 
which is shown as black square at the bottom of Fig. [1] 

The same QCP in Ref. l4ll (note factor of i/2 instead of 
t in their Hamiltonian) was estimated to be {U/t)c ^ 7, 
substantially different from the result we obtain here^. 
However, for p = 0.75 Ref. [lO (note factor of t/2 instead 
of t in their Hamiltonian) finds {U/t)c ^ 10 and in the 
hard-core limit Ref. 44 (note factor of U in stead of U/t in 
their Hamiltonian) finds {U/t)c = 9.9 ± 0.4 with p = 0.5 
in reasonable agreement with our results. We note that 
in Ref. 41 the location of the QCP was obtained partly by 
requiring the frequency dependent conductivity, (t{uj), to 
scale with w at the critical point in violation of Eq. ([3|) . 
We now turn to a discussion of our results for the dynamic 
conductivity. 

As discussed above, the dynamic conductivity calcu- 
lated with Matsubara frequency is expected to satisfy 
the scaling relation Eq. ([3]). The implication of this scal- 
ing law for numerical simulations is that cr(a;, T, L) in the 
limit L — >■ oo as T approaches should be a function of a 
single variable huj / (ksT). Therefore, if the dynamic con- 
ductivity curves calculated from different lattice sizes are 
plotted versus w, one would expect deviations from a sin- 
gle curve to be visible at the smallest Matsubara frequen- 
cies since finite size dependence there is likely the smallest 
and the data should be close to scaling with hoj/jkBT) 
even without extrapolation to L —>■ oo. See also Ref. l39l . 
Such behavior is clearly visible in Fig. [3 where SSE re- 
sults for the dynamic conductivity a{ujk)/crQ at the crit- 
ical point Zt/Uc = 0.182 are shown. All curves from 
4 different lattice sizes overlap with each other at high 
frequency limit but differ at the low frequency side. A 
rough estimate of the high frequency universal conduc- 
tivity, a* using the lattice size L = 16 gives a* ~ 0.17o'q. 
This is in surprisingly good agreement with simulations 



S, 




FIG. 7: (Color online) Dynamic conductivity scaling plot 
a{ujk)/o-Q vs Matsubara frequency Uk for Zt/Uc ~ 0.182 and 
i = 0.5. 



of the quantum rotor model, Eq. ([3]), where one finds 
a* — (0.14 ± 0.01)tTQ, — as weU as with the exact di- 



8 



agonalization result (0.15 ± 0.01)(tq for a hard-core BH 
modelj^ as one would expect for a universal quantity. 
However, this value is much less than the value found by 
Batrouni et al. (0.47±0.08)ctq.'''^ The discrepancy could 
be due to an inaccurate estimate of the QCP in Ref. ^ 
obtained partly by assuming that (t(w) should scale with 
Lu at the quantum critical point. 
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FIG. 8: (Color online), (a) Scaling curves of the superfluid 
density from grand canonical simulations at fi/U = 0.375. 
The inset shows the slope of these curves at the QCP on a 
log-log scale, (b) The compressibility k in the same region. 



In addition we have performed canonical simulations 
with p = 1, and determined the QCP at [(Zi/[/)c = 
0.108(2)], which is slightly different from the fixed chem- 
ical potential (/x/C/ = 0.375) analysis value of {Zt/U)c — 
0.112(1) to be discussed in the next section. 



B. Grand canonical simulations ^/U — 0,0.375, 1 

Close to the commensurate filling of p = 1 correspond- 
ing to u/U = 0.375 the QCP has been determined in 
Ref. [3^ for the clean BH model (see also solid triangle in 
Fig. [1]). Here we are interested in further investigating 
the effect of disorder on this transition by keeping all the 



other parameters the same as in Ref. |39| while setting 
the disorder potential amplitude A = U. 

Note that in generating the disordered onsite potential 
Ci for the system, we make sure that the average of the 
disordered onsite potential ^ ■ ei/N ~ for each disor- 
dered configuration of the system. This can be achieved 
by first generating a disorder configuration and calcu- 
lating the average of the disordered potential, then sub- 
tracting this average from the generated potential at each 
lattice site. 

Our results are shown in Fig. [SJa) where L^[< ps >] is 
plotted for different system sizes close to the QCP. The 
scaling law here is the same as outlined above for the 
simulations in the canonical ensemble and again we note 
that the relation z = d was assumed. As can be seen 
from Fig. [UJa) the curves cross very nearly in a single 
point and it is easy to estimate the corresponding QCP 
points as {Zt/U)c = 0.112(1). Very interestingly, when 
this is compared to the results for the clean system QCP 
of {Zt/U)c = 0.2385 from Ref. \M, it is clear that the 
introduction of disorder at this commensurate filling has 
enhanced the superfluid region. The QCP with (■) and 
without (a) disorder at fj,/U = 0.375 are shown along the 
dashed line in Fig. [T] A similar indication for disorder- 
enhanced superfluidity has been found by Krauth et al. 
fi^ but the effect is even more noticeable in our results. 

Also shown in Fig.[51[b) (as an inset) is the slope of L'^Ps 
at the QCP. From the scaling relation Eq. (|14p it is clear 
that this slope should scale with L^/"^. The solid line in 
the inset is a flt to the data yielding t^ ^ 1 in agreement 
with the quantum version of the Harris criterion v > 2/d. 

In Fig. ^h) are shown results for the compressibility, 
K, through the QCP. As can be seen, the compressibility 
remains constant through the QCP. It is expected that 
in the presence of disorder the compressibility will scale 
as K ~ §'^{d,-z)^ rpj-^g ^^^^ ihaX our results for k remain 
constant through the transition is consistent with the re- 
lation z — d, although it does not rule out z < d. 

At the QCP we also measure the dynamic conductiv- 
ity with grand canonical QMC simulations, shown in Fig. 
[9l Like the simulations in the canonical ensemble, we 
observe that the data appear to scale with ujk at high 
frequency but very clearly deviate from this behavior at 
low frequencies. Again we emphasize that this is con- 
sistent with a{uj) scaling with huj/{kBT) rather than uj. 
Here we can also attempt to estimate a* and we flnd it to 
be around 0.17(7q, in very good agreement with the re- 
sults for the canonical simulations from the previous sec- 
tion as well as with the simulations of the quantum rotor 
modelJ^ As before, we note that the universal conduc- 
tivity, a* estimated here and in the previous section cor- 
responds to the conductivity value in the high frequency 
limit, not the DC conductivity measured experimentally. 

We have also performed simulations ai p/U = 1 where 
we flnd {Zt/U)c = 0.101 as well as for p/U — where 
we find {Zt/U)c = 0.192. These QCPs are also shown as 
solid black squares in Fig. [T] 
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FIG. 9: (Color online) Dynamic conductivity a{u)k)/oQ vs 
Matsubara frequency ujk for Ztc/U = 0.112 in a disordered 
grand canonical system at ^l/U — 0.375. 



IV. CONCLUSION 



results we present here are direct simulations of the full 
BH model with large disorder. Our results clearly indi- 
cate that these three models in the presence of disorder 
are in the same universality class. 

At the QCP, we also compute the dynamic conductiv- 
ity as a function of the Matsubara frequency. The uni- 
versal conductivity at high frequency limit is estimated 
to be around O.ITctq in good agreement with previous 
calculations on quantum rotor models.— Most notably, 
at the QCP, the dynamic conductivity shows clear devi- 
ations from scaling with w, consistent instead with the 
expected scaling with w/(fcsT). In addition, we checked 
in detail the convergence with the number of disorder 
reahzations finding that a large number of disorder re- 
alizations (usually on the order of 10*) are necessary to 
have a fully convergent superfluid density value. This 
implies that some previous calculations in the literature 
with disorder realizations of the order of 10^, might not 
be reliable. This has allowed us to determine the QCP for 
several values of fi/U for the disordered Bose-Hubbard 
model. 
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